clear
clc
close all

%初始已知控制点
x_control=0:0.1:3;
y_congrol=sin(2*pi*0.5*x_control)+cos(2*pi*0.8*x_control)+2+0.02*rand(size(x_control));

x=(-1:0.01:4);
y=RBF1(x_control,y_congrol,x,'gaussian',0.3,0,0.001);

figure()
hold on
plot(x_control,y_congrol,'o')
plot(x,y)
hold off

function y=RBF1(x0,y0,x,Method,Rs,Npoly,Error)
%一维RBF插值，输入x0散点，y0值，是已知控制点。x是插值点。
%Method方法，默认'linear'。
%'linear',|R|
%'gaussian',exp(-(R/Rs)^2)
%'thin_plate',R^2*log(R)
%'linearEpsR',|R|。分段线性插值，要求s更大一些
%'cubic',|R^3|

%Rs，插值核作用半径。Rs对于'linear'，'cubic'，'thin_plate'无影响。Rs大概和点和点之间的距离差不多就行。
%Npoly多项式拟合。默认是1，只拟合1次项。
%Error误差。在(-∞,∞)区间。默认是0，无误差。表示可以在部分误差范围内去插值。
%Error一般在0~1之内。如果只是为了收敛，可以比较小，在0.1以内就可以有很好效果；如果想实现平滑，可适当增大，大于1。

%示例：x0=0:0.3:pi;y0=sin(2*x0)+0.2*x0;x=-1:0.01:4;y=RBF1(x0,y0,x,'linear',0.2,1,0.0);
%示例：x0=0:0.1:3;y0=0.25*x0.^2-0.5*x0+0.1*rand(size(x0));x=-1:0.01:4;y=RBF1(x0,y0,x,'gaussian',0.3,2,0.01);

N=length(x0);%点的数量，也是RBF基函数的数量
%整理为列向量
x0=x0(:);
y0=y0(:);
x=x(:);

%函数默认信息
if nargin<4 || isempty(Method)
    Method='linear';%默认线性核函数
elseif nargin<5 || isempty(Rs)
    Rs=1.1*(max(x0)-min(x0))/N;%假设空间均匀分布
elseif nargin<6 || isempty(Npoly)
    Npoly=1;%默认拟合1次项
elseif nargin<7 || isempty(Error)
    Error=0;%默认输入数据没有误差
end

%选择核函数
K_method=0;
switch Method
    case 'linear'
        K_method=1;
        fun=@(RMat) Kernel_Linear(RMat,Rs);
    case 'gaussian'
        K_method=2;
        fun=@(RMat) Kernel_Gaussian(RMat,Rs);
        Error=-Error;%因为零点不是零，远点是0，所以这里误差为负
    case 'thin_plate'
        K_method=3;
        fun=@(RMat) Kernel_Thin_plate(RMat,Rs);
    case 'linearEpsR'
        K_method=4;
        fun=@(RMat) Kernel_LinearEpsR(RMat,Rs);
        Error=-Error;%因为零点不是零，远点是0，所以这里误差为负
end

%将原始数据分离出多项式项Npoly
if Npoly==0
    C=ones(N,1);%常数项
    wC=mean(y0);
elseif Npoly==1
    C=[ones(N,1),x0];%常数项+一次项
    wC=C\y0;
elseif Npoly==2
    C=[ones(N,1),x0,x0.^2];%常数项+一次项+二次项
    wC=C\y0;
else
    error('只支持Npoly=0,1,2')
end
yC=C*wC;%多项式项
y1=y0-yC;

%计算每个基函数在各个点的值
K2=zeros(N);%每一列对应一个函数
%计算距离矩阵
DisMat=squareform(pdist(x0));
K3=fun(DisMat);%每一个核函数的中心点在节点上
K3=K3-Error*eye(N);%把误差函数加入

w=K3\y1;%计算权重

%根据权重计算插值
y=zeros(size(x));
for k=1:N %计算每个核的贡献，然后叠加
    R_k=abs(x-x0(k));
    yt=w(k)*fun(R_k );
    %plot(t,yt);
    y=y+yt;
end

NOut=length(y);
%再还原回多项式项
if Npoly==0
    C=ones(NOut,1);%常数项
elseif Npoly==1
    C=[ones(NOut,1),x];%常数项+一次项
elseif Npoly==2
    C=[ones(NOut,1),x,x.^2];%常数项+一次项+二次项
end
y=y+C*wC;%再加上多项式项

%检查结果
if max(abs(w))/max(abs(y1))>1e3
    warning('结果未收敛，建议调整误差Error');
elseif (max(y)-min(y))>5*(max(y0)-min(y0))
    warning('结果未收敛，建议增大区间Eps，或者调整误差Error');
end

end

function z=Kernel_Linear(R,s)
%R距离
z=abs(R);%线性函数
end

function z=Kernel_Gaussian(R,s)
%R距离
%s大概和采样点间距差不多就行，可以略大（更胖）
z=exp(-R.^2/2/s^2);%正态函数
end

function z=Kernel_Thin_plate(R,s)
%R距离
z=R.^2.*log(R);%thin_plate
end

function z=Kernel_LinearEpsR(R,s)
%R距离
%s大概和采样点间距差不多就行，可以略大（更胖）
s=2*s;%这里要更大
z=s-R;%线性函数
z(z<0)=0;
end
